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INTRODUCTION 

Due  to  recent  economic  developments  such  as  the  energy 
crisis  and  persistent  cost  push  inflation  there  has  been  consid- 
erable interest  or  the  structure  of  national  economies,  and  how 
crisis  effects  economic  activities.  One  method  that  is  used  to 
study  economic  systems  is  Input  Output  (10)  analysis  which  is 
particularly  powerful  in  the  study  of  interindustry  activity. 
The  purpose  of  this  paper  is  to  discuss  the  solution  of  the  10 
equation  using  several  methods,  the  goal  being  to  develop  a 
method  to  solve  a  nonlinear  Input  Output  model  which  will  partly 
alleviate  some  of  the  restrictions  of  linear  Input  Output 
analysis. 

Input  Output  analysis  is  an  econometric  method  which  at- 
tempts to  explain  all  industrial  activity  by  a  simple  cause  ef- 
fect relationship.  The  first  critical  assumption  of  10  analysis 
is  that  all  goods  in  a  product  group  are  manufactured  in  an 
identical  manner.  From  this  point  on  when  the  term  good  is  men- 
tioned, it  is  to  be  taken  as  some  representative  good  in  one  of 
the  product  groups.  The  amount  of  a  good  that  society  needs  to 
produce   is   the  amount  to  be  supplied  for  final  demand   plus  the 


By  final  demand  we  mean  several  things:  goods 
consumed  by  households  and  government,  goods  sold  for 
export,  and  also  goods  used  for  investment.  Most  logically 
investment  would  be  treated  as  an  input  to  production,  but 
investment  can  be  a  very  nonlinear  function  of  output.  Thus 
in  general  economists  find  it  easier  to  determine  investment 
demands  exogeneously .  It  is  possible  that  the 
representation  ot  investment  as  a  nonlinear  function  of 
output  could  make  the  endogenous  determination  of  investment 
demands  more  realistic. 


amount  used  in  the  production  of  other  goods.  The  second  criti- 
cal assumption  of  10  analysis  is  the  following:  when  a  good  is 
used  in  the  production  of  another  good,  the  amount  used  in  this 
activity  is  always  in  a  fixed  proportion  to  the  total  production 
of  the  other  good.  These  proportionality  constants  are  termed 
the  technical  coefficients,  and  since  in  general  it  is  possible 
that  all  goods  can  be  used  directly  in  the  production  of  all  oth- 
er  goods,  if  we  have  an  n  good  economy  then  there  are  n*  techni- 
cal coefficients.  Vvhen  the  technical  coefficients  are  arranged 
in  a  nxn  matrix,  this  matrix  is  called  the  technical  coefficients 
matrix  or  A  matrix,  though  this  should  not  be  confused  with  the  A 
matrix  found  in  control  systems  literature.  If  our  unit  of  value 
is  dollars  then  the  element  a.  .  represent  the  dollar  value  of  the 
i  good  required  in  the  direct  production  of  one  dollar  of  the 
j    good.   Thus  the  a^'s  are  positive  fractions.   The  sum  of  the 

t"  h 

j  column  of  the  A  matrix  represents  the  fraction  of  the  total 
cost  of  producing  the  good  j  that  is  embodied  in  the  goods  used 
in  direct  production  of  the  j    good.   Then 


v,  =  1  -  I     a-,  (0.1) 

J       i=i   *J 


represents  the  unit  cost  of  production  not  embodied  in  the  goods 
used  in  the  direct  production  of  the  j  good.  V.  is  termed  the 
value  added  for  good  j  and  it  consists  of  labor  costs,  interest, 
(on  the  capital  used  in  the  direct  production  of  the  jtn  good,) 
direct  business  taxes,  and  profits. 


If  y^  is  the  amount  of  the  i  good  sold  to  final  demand, 
then  the  equation  which  represents  the  total  production  of  the 
i   good,  that  is,  x.,  is 


*i  s  yt  «■   2  a..x,  .  (0.2) 


The  summation  term  in  (0.2)  is  the  amount  of  good  i  needed  in 
the  direct  production  of  all  goods.  Proceeding  in  the  same  way 
for  the  other  n-1  goods,  the  total  demand  for  all  goods  is  the 
solution  of  the  matrix  equation 

Ax  +  y  =  x   or   (I  -  A)x  =  y  .  (0.3) 

We  reiterate  the  basic  assumptions  of  10  analysis;  1)  that  all 
goods  in  the  same  product  class  are  assumed  to  be  made  in  the 
same  way,  2)  that  the  amount  of  an  input  good  used  in  the  direct 
production  of  another  good  is  always  in  a  fixed  proportion. 

The  first  question  that  needs  to  be   answered   is  whether 
solutions   to   (0.3)   do   indeed  exist.   But  since  we  also  desire 

that  solutions  must  correspond  to  actual  economic   behavior,   the 

2 
solution   vectors  should  be  positive   if  the  final  demand  vector 

is  positive.  Bellman [1,  pp.296,  Theorem  6]  has  shown  that  if   the 

column  sums  of  A  are  less  that  1,  then  there  is  a  unique  positive 


o 

A  vector  is   said   to   be   positive   if   all   of   its 

elements  are  positive. 


solution  vector  for  each  positive  right  hand  side.    Furthermore, 
the  eigenvalues  of  A  lie  inside  the  unit  circle,  and 


(I  -  A)"1  =  I  +   2  A1  ,  (0.4) 


i  =  l 


with  lim  Am  =  0,  e.g.  see  Isaacson  and  Keller [15,  pp.15,   Theorem 

nHro 
5].    If   k   terms  are  used  to  approximate  (I  -  A)"1  then  (k-l)n 

multiplications  must  be  performed. 

In  the  past,  eg.  Chenary  and  Clark[17],  (I-A)  was  crude- 
ly approximated  by  this  method.  The  objective  of  this  paper  is 
to  discuss  a  method  by  which  the  10  equation  can  be  solved  with 
far  more  efficiency  and  accuracy.  By  a  factorization  method 
which  will  be  discussed  in  Chapter  1,  only  J:n^  multiplications 
are  required  to  completely  factorize  the  matrix  into  a  product  of 

triangular  matrices,  which  can  then  be  solved  for  arbitrary  right 

2 
hand  sides  in  n  multiplications.   The  first  section  of  Chapter  1 

will  also  analyze  the  numerical  stability  properties  of  factori- 
zation of  Input  Output  matrices  by  elementary  transformations, 
(sometimes  denoted  as  LU  decomposition,)  and  will  reexamine  the 
property  of  the  A  matrix  so  that  (0.3)  will  admit  only  positive 
solutions  for  arbitrary  positive  final  demand  vectors.  The 
second   section   in   this   chapter  proposes  a  simple  method  which 
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It  is  standard  practice  in  numerical  analysis  to  only 

count  the  number  of  multiplications  or  divisions  in  a 
computation  since  they  are  more  time  consuming  than  addition 
or  subtraction  operations. 


will  allow  modifications  of  particular  rows  and  columns  of  the  A 
matrix  that  will  not  require  the  complete  ref actor ization  of  the 
matrix.  This  method  has  a  very  useful  property  that  allows  the 
successive  updating  of  the  10  matrix  without  any  algorithmic  com- 
plexity. Procedures  of  this  type  are  generally  referred  to  as 
factorization  modification,  and  several  papers  have  been  written 
on  this  subject,  including  Gill  and  Murray[13],  Golub  et  al.[12], 
Sameh  and  Bezdek[2],  and  Noh  and  Sameh[3].  But  with  the  excep- 
tion of  Gill  and  Murray[13],  these  only  deal  with  factorization 
by  orthogonal  transformations,  which  will  be  more  numerically 
stable  for  general  matrices,  but  require  more  storage  and  compu- 
tation than  factorization  by  elementary  transformations  thus  mak- 
ing them  comparably  more  expensive.  By  examining  the  structure 
of  10  matrices,  we  shall  see  that  the  use  of  elementary  transfor- 
mations in  the  factorization  of  10  matrices  is  quite  appropriate. 
We  do  not  mean  to  distract  the  reader  away  from  orthogonal 
transformation  factorization  methods.  As  Bierman[5]  points  out, 
(this  article  is  an  excellent  survey  of  numerical  techniques  in 
control  and  other  applications,)  in  general  these  methods  lend  to 
many  algorithmic  advantages,  and  the  numerical  stability  which 
results  is  important  in  control  problems  which  can  be  ill- 
conditioned. 

The  next  section  discusses  the  updating  of  solutions  of 
linear  equations  when  only  a  few  elements  of  several  columns  or 
rows  are  changed.  We  will  denote  this  method  by  "solution  per- 
turbation" since  the  method  depends  on  solutions  to  a  nominal 
system  of  equations.   This  is  a  bit  of  a  misnomer,  since   pertur- 


bation  implies  small  changes.  However  here  the  changes  in  ele- 
ments of  the  matrix  need  not  be  small  in  any  way.  In  many  cases 
this  method  has  a  tremendous  computational  advantage  over  factor- 
ization modification  but  it  does  have  its  drawbacks.  This  method 
requires  that  columns  of  the  nominal  inverse,  (that  is  of 
(I  -  A)  before  any  elements  have  changed,)  to  be  computed.  In 
the  case  of  modifying  an  entire  column  solution  perturbation 
would  require  that  the  entire  nominal  inverse  be  computed.  Since 
the  computation  of  an  inverse  requires  approximately  three  times 
the  computation  required  to  factorize  a  matrix,  this  method  is 
not  appropriate  if  many  elements  in  a  column  are  to  be  modified. 
Also  the  method  has  a  disadvantage  in  that  several  successive  up- 
dates of  the  solution  due  to  new  perturbations  in  the  matrix  ele- 
ments are  difficult  to  perform  since  there  is  no  efficient  way  to 
compute  the  updated  inverses. 

In  Chapter  3  a  nonlinear  10  model  is  proposed  which  can  be 
solved  quite  efficiently  with  the  eclectic  use  the  algorithms 
described  in  the  preceding  chapters.  This  approach  largely  elim- 
inates the  linearity  contraint  of  the  standard  10  model  while  the 
extra  cost  of  solving  the  nonlinear  model  is  quite  small.  Using 
solution  perturbation  the  nonlinear  equation  that  must  be  solved 
iteratively  is  reduced  to  smaller  order.  Once  the  residuals  are 
sufficiently  small  factorization  modification  will  be  utilized  to 
find  the  complete  solution.  These  algorithms  have  been  coded  on 
a  minicomputer  system,  and  the  computation  of  solutions  of 
350  -order  10  equations  is  quite  feasible  on  such  a  computer 
system  when  these  algorithms  are  used. 


CHAPTER  1 
FACTORIZATION  MODIFICATIONS 
SECTION  1.1 

MATRIX  DECOMPOSITION  BY  ELEMENTARY  TRANSFORMATIONS,  AND 
SIMPLIFICATIONS  THEREOF  ON  10  MATRICES. 


This  section  describes  the  LU  decomposition  algorithm. 
The  presentation  will  also  include  a  discussion  of  properties  of 
10  matrices  that  will  simplify  the  algorithms  to  be  described  in 
the  next  section. 

As  with  all  decomposition  algorithms,  in  the  decomposition 
of  the  matrix  B,  (in  our  case  B=  I-  A,  where  I  is  the  identity 
matrix,  and  A  is  the  technical  coefficient  matrix  defined  in  the 
introductory  section,  thus  b^.  <  0  and  bii  =  1  -  aii  >  0,)  a 
transformation  0  is  determined  such  that 

QB=  U  (1.1.1) 

where  U  is  a  upper  triangular  matrix.  To  solve  the  system  of 
equations 

Bx  =  y  (1.1.2) 

we  premulitply  both  sides  of  (1.1.2)  by  Q,  so  that 

QBx  =  Ux  =  Qy  (1.1.3) 

and  the  equation  Ux  =  Qy  can  be  solved  by  back  substitution.  In 
the  pure  form  of  LU  decomposition  Q  is  the  composition  of  only 
elementary  transformations.   An  elementary  transformation   has   a 


8 

simple  matrix  representation.  Its  diagonal  elements  are  all  one 
and  it  has  only  one  nonzero  off  diagonal  element  which  we  shall 
call  the  multiplier  element.  If  Mi-  is  an  elementary  transforma- 
tion then  it  looks  like: 


m 


ID 


(1.1-4) 

To  simplify  the  notation,  the  subscripts  below  a  symbol 
representing  an  elementary  transformation  will  denote  the  posi- 
tion of  the  multiplier  element  in  the  elementary  transformation 
matrix.  Thus  m.  .  is  in  the  position  (i,j)  of  the  matrix.  It  is 
easily  verified  the  the  inverse  of  an  elementary  transformation 
is  merely  the  elementary  transformation  with  the  nonzero  off  di- 
agonal element  of  opposite  sign.  The  direct  multiplication  of 
elementary  matrices, 


N.  =  to      .m_  t  _: M^,0  -iM-,,  • 


(1.1.5) 


looks  like: 


j  +  2,j 


m 


no 


(1.1.6) 


The  LU  decomposition  algorithm  determines  the  transforma- 


tion 


Q  =  Nn-lNn-2 


N2N1 


(1.1.7) 


If  for  each  j  we  define 


Nj(Ej)=  Nj(Nj.1Nj_2...N1B)  , 


(1.1.8) 


then  the  algorithm  determines  N.  such  that  it  zeroes  the  elements 
of  the  j   column  of  B.  below  the  diagonal  element.   Thus 


m 


13 


13 


33 


for  j  +  Ki<n 


(1.1.9) 
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where  the  n^j's  are  elements  of  (1.1.6),  B^  =  (E^),  and  d".^  is 
called  a  pivot  element. 

A   reader   familiar   with   decomposition    by   elementary 

4 
transformations   probably  notes  that  row  permutation  has  been  om- 
itted.  The  reason  for  this  will  be  clear  shortly. 

Proposition ;  If  all  the  column  sums  of  A  are  less  than  unity  then 
the  pivot  elements  are  nonzero. 

Proof :  Clearly  if  the  column  sums  of  A  are  less  than  unity  then 
the  same  must  hold  for  the  i  leading  principal  minor  of  A. 
Thus  by  Bellman [1,  pp.296,  theorem  6]  the  principal  minors  of 
(I-A)  are  nonsingular.  Stewart[10,  pp.120,  theorem  2.5]  has 
shown  that  if  the  principal  minors  of  a  matrix  are  nonzero  then 
if  the  matrix  is  decomposed  by  elementary  transformations  without 
pivoting,  the  pivot  elements  will  be  nonzero. 

There  is  a  fundamental  result  in  Input  Output  analysis 
which  describes  the  necessary  and  sufficient  condition  on  A  so 
that  the  10  equation  admits  only  positive  solutions  vectors  for 
positive  final  demand  vectors.  This  is  called  the  Simon-Hawkins 
condition [9] .  The  condition  is  that  all  the  principal  minors  of 
(I-A)  must  be  positive.  Let  us  examine  the  computation  of  the 
solution  of  10  equation  by  factorization  so  that  we  can  see  when 


4 

See   Forsythe   and   Moler[16],    or    Isaacson    and 

Keller [15]    on    introductory   material   on   decomposition 

methods. 
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it  is  impossible  for  a  solution  to  be  a  nonpositive  vector  when 
the  final  demand  vector  is  positive.  After  applying  the 
transformation  Q  to  the  final  demand  vector,  since  all  the  multi- 
plier elements  of  Q  are  nonnegative,  each  element  of  Qy  must  be 
greater  than  or  equal  to  the  corresponding  element  of  y.  In  the 
solution  of  Ux=Qy,  since  the  off  diagonal  elements  are  all  nonpo- 
sitive, (thus  in  the  back  substitution  all  sums  are  on  numbers  of 
nonnegative  sign,)  the  only  way  that  a  negative  element  can  occur 
in  the  solution  is  if  a  diagonal  element  of  U  is  negative.  Thus 
if  all  the  diagonal  elements  of  U  are  positive,  (which  implies 
that  all  of  the  multiplier  elements  of  the  elementary  transforma- 
tions must  be  positive,)  then  it  is  impossible  for  the  factorized 
equation  to  admit  a  nonpositive  solution  vector  for  a  positive 
final  demand  vector.  This  is  exactly  how  the  Simon-Hawkins  con- 
dition is  checked:  the  determinant  of  the  k  principal  minor  of 
(I  -  A)  is  merely  the  product  of  the  k  topmost  elements  of  the 
diagonal  of  U  since  det|Q|=l. 

In  the  definition  of  the  transformation  Q  one  order  of  ap- 
plying the  elementary  transformations  was  given.  Since  elementa- 
ry transformations  are  nearly  identity  matrices  one  would  expect 
that  these  transformations  will  commute  in  certain  cases,  (clear- 
ly any  permutation  of  the  elementary  transformations  in  N.  will 
yield  the  same  transformation.)  The  elementary  transformations 
can  be  commuted  just  as  long  as  no  element  becomes  nonzero  that 
was  intentionally  zeroed  by  the  application  of  a  previous  elemen- 
tary transformation.   Thus: 


12 


Q  *  *1nn-l---Mn2"-M32Mnl---M21 


(1.1.10) 


=   Mnn-1"-M43^42K41M32M31M21 


(1.1.11) 


If  we  apply  Q  by  (1.1.10)  then  after  the  first  n-1  elementary 
transformations  are  applied,  the  partially  decomposed  matrix  has 
the  structure: 


X  X  X  X  X  X  X 

0  X  X  X  X  X  X 

0  X  X  X  X  X  X 


0 


0 


u 


(1.1.12) 


After    the  application   of   the  next   n-2   elementary      transformations 
the  matrix   has   the   structure: 


0  0 
0  0 
0        0 


XXX 


0  X  X  X  X  X  X 

0  0  X  X  X  X  X 

0  0  X  X  X  X  X 


X  X  X  X  X 

X  X  X  X  X 


(1.1.13) 


If  we  apply  Q  by  (1.1.11)  then  after  the  application  of  the  first 
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elementary  transformation  the  partially  decomposed  matrix  has  the 
structure : 


x   x   x    x    x   x   x 

0     X     X     X     X     X     X 


(1.1.14) 


Then  after    the  application   of   the  next   two  elementary   transforma- 
tions  the   matrix   has   the    structure: 


X  X  X  X  X  X  X 
0  X  X  X  X  X  X 
0     0     X     X     X     X     X 


(1.1.15) 


Another  ordering  which  will  reduce  page  faults  in  virtual 
computers  or  10  requests  in  successively  reading  in  and  writing 
out  parts  of  the  matrix  on  machines  with  little  main  memory  in- 
volves storing  the  matrix  in  blocks  of  rows.  The  topmost  par- 
tially decomposed  block  is  completely  decomposed.  Then  using  this 
completely  decomposed   block   the  remaining  blocks  are  partially 


14 

decomposed  using  this  block,  that  is,  all  the  transformations 
that  need  this  completely  decomposed  block  will  be  applied  to  the 
lower  blocks.  Thus  this  block  will  no  longer  be  needed.  This 
will  reduce  page  faults  or  10  requests  by  a  factor  equal  to  the 
number  of  rows  that  are  stored  in  each  block.  All  these  methods 
are  equivalent  not  only  mathematically  but  numerically,  that  is, 
if  one  is  careful  about  the  ordering  of  the  computations,  all  of 
the  methods  will  achieve  identical  matrices  in  terms  of  the  bit 
patterns. 

The  name  of  the  method,  LU  decomposition,   refers   to   the 
definition 

LU  =  B  (1.1.16) 

where  L  is  lower  triangular.  By  inspection 


L  =  Q"1  .  (1.1.17) 


This  is  easily  shown  by  the  direct  multiplication  of  the  inverses 

of  the  elementary  transformations. 

3 

The  algorithm  performs  approximately  n  /3  multiplications 

and  requires  only  the  original  matrix  as  work  space.  The  algo- 
rithm is  the  fastest  and  most  compact  of  all  decomposition  algo- 
rithms, Moler [ 11] . 
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SECTION  1.2 

MULTIPLE  ROW  AND  COLUMN  MODIFICATIONS 

In  this  section  an  algorithm  will  pp  discussed  which  will 
allow  the  modification  of  the  factorization  of  a  matrix  when  cer- 
tain preselected  rows  and  columns  of  the  matrix  are  changed.  A 
significant  advantage  of  this  algorithm  is  that  each  modification 
of  the  factorization  adds  no  complexity  to  solving  the  new  system 
of  equations,  since  the  modified  factorization  is  solved  in  the 
same  manner  as  the  original  system.  Also  the  modified  factoriza- 
tion is  identical  to  the  factorization  resulting  from  completely 
factorizing  the  modified  matrix.  As  we  shall  see  the  computation 
of  the  modified  factorization  and  the  computation  of  solutions  to 
the  modified  system  of  equations  can  be  performed  simultaneously. 
This  is  useful  when  the  computations  are  performed  on  a  minicom- 
puter. Experience  has  shown  that  the  cost  of  reading  in  the  ma- 
trices is  the  most  significant  cost  in  the  solution  of  a  large 
factorized  system  of  equations.  Thus  streamlining  the  algorithms 
with  respect  to  the  number  of  times  that  the  matrices  must  be 
read  in  from  secondary  storage  is  worthwhile.  In  Aoki[4]  and 
Householder [5] ,  LU  decomposition  is  introduced  in  a  manner  that 
would  lend  itself  to  factorization  modification  methods,  though 
there  was  no  intention  to  discuss  the  subject. 

The  algorithm  for  factorization  modification  greatly  sim- 
plifies if  the  rows  to  be  modified  are  at  the  bottom  of  the  ma- 
trix, and  the  columns  are  on  the  right  hand  border.  Thus  before 
performing  the  factorization  we  exchange  the  rows  to  be  modified 
with  the  bottom  rows  of  the  matrix,  and  exchange  the   columns  to 
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be  modified  with  the  columns  on  the  right  hand  border.  Let  P, 
P~  «P  ,  be  the  transformation  that  exchanges  the  rows  and 
columns.  Clearly  the  matrix  representation  of  P  has  n^-n  zero 
elements,  the  other  n  elements  being  unity.  Rows  of  the 
transformation  corresponding  to  coordinates  that  are  not  permuted 
have  1  on  the  diagonal.  If  coordinate  i  is  to  be  exchanged  with 
the  j  coordinate,  then  the  elements  in  positions  (i,j)  and 
(j,i)  would  be  one.  Therefore  for  all  i,j  p.  .«p..  or  P*Pt.  Also 
it  is  easily  verified  that  PP=I. 

For  example  let  there  be  two  industries  for  which  we 
desire  column  or  row  modifications,  their  positions  being  say  1 
and  3.  We  wish  to  construct  a  transformation  with  the  properties 
above  which  permutes  these  industries  so  that  their  rows  are  at 
the  bottom  and  the  columns  are  on  the  right  hand  border.  P  would 
then  look  like: 


0 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

(1.2.1) 


The  above  example  verifies  that  P   is   syntpetric.    The   permuted 
form  of  B  will  be  denoted  by 


B  ■  P  BP  -  PBP  . 


(1.2.2) 
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Let  Q  be  the  composition  of  M.,'s  for  l<j<n-k,  j+l<i<n. 
Q  will  partially  decompose  the  matrix  B.  The  elementary 
transformations  are  to  applied  as  not  to  destroy  zeros  previously 
introduced.   The  partial  decomposition  will  be  denoted  by 


QpB  *  0  . 


(1.2.3) 


In  our  6x6  example  above,  pre-  and  post-multiplying  the  original 
matrix  A  by  P  will  permute  the  rows  and  columns  1  and  3  to  the 
bottom  and  right  border  of  the  matrix.  Note  that  P(I  -  A) P  ■ 
(I  -  PAP) .  After  Q  is  applied  to  our  permuted  (I-A)  matrix,  the 
resultant  matrix  0  has  the  structure: 


0 

0 

X 

X 

X 

X 

0 

0 

0 

X 

X 

X 

0 

0 

0 

0 

X 

X 

0 

0 

0 

0 

X 

X 

(1.2.4) 


At  this  point  any  of  the  last  k  columns  of  B  can  be  modi- 
fied. Let  us  assume  that  we  want  to  replace  a  column  of  B  which 
corresponds  to  one  of  the  last  k  columns  of  B  by  a  column  vector 
g.  0  is  updated  merely  by  replacing  the  corresponding  column  of 
u*  by  QpPg.  Changing  more  such  columns  requires  the  identical  pro- 
cedure  for   each   column.   The  transformation  of  g  and  the  right 
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hand  side  vector  by  Q  p  can  be  performed  simultaneously.  This 
feature  in  very  useful  when  the  computation  is  performed  on  a 
minicomputer  system. 

Also  any  of  the  last  k  rows  of  B  can  be  modified.  If  a 
row  of  B  which  corresponds  to  one  of  the  last  k  rows  of  B  is  to 
be  replaced  by  the  row  vector  h  ,  then  the  corresponding  elemen- 
tary transformations  in  Q  which  zeroed  out  the  first  n-k  ele- 
ments  of  corresponding  row  of  0  are  recomputed.  Suppose  that  in 
our  6x6  example  we  desire  to  replace  the  first  row  of  B  by  h  . 
This  would  correspond  to  replacing  the  5  row  of  B  by  (Ph)  . 
Thus  we  recompute  the  elementary  transformations  M_,  through  Mcd. 
Before  these  elementary  transformations  are  applied,  the  modified 
0  has  the  structure: 


r 


X 

X 

X 

X 

X 

X 

0 

X 

X 

X 

X 

X 

0 

0 

X 

X 

X 

X 

0 

0 

0 

X 

X 

X 

X     X     X     X     X     X 

0     0     0     0     X     X 


(1.2.5) 


Note  that  this  procedure  does  not  really  violate  the  rule  that  no 
element  that  has  been  intentionally  zeroed  be  set  nonzero  by  a 
elementary  transformation  applied  later.  Since  we  are  recomputing 
all  the  elementary  transformations  that  zeroed  out  the  elements 
in  the  row  that  is  to  be  modified,  it  is  as  though  these  elemen- 
tary  transformations  were  never  applied  before  the  modification 


19 

was  made. 

For  any  number  ot  row  and/or  column  modifications,  the 
modified  factorization  will  be  identical  to  the  factorized  matrix 
which  would  result  from  partially  decomposing  B,  where  B  is  the 
matrix  with  the  row  and/or  column  modifications.  In  fact  when 
the  original  rows  and  columns  are  replaced  in  the  factorization, 
the  new  file  is  exactly  bit  comparable  to  the  original.  To 
achieve  this  the  reader  is  warned  that  care  must  be  taken  in  the 
ordering  of  operations  in  the  original  decomposition  and  the  rows 
and  columns  modification  algorithms. 

To  solve  the  original  or  a  modified  system  of  equations, 
the  rest  of  the  elementary  transformations  are  applied  that  will 
reduce  0  to  a  triangular  matrix.  This  second  composition  of  ele- 
mentary transformations,  (the  M^s'  where  n-k+l<i<n-l  and 
i+l<j<n,)  are  applied  in  an  order  that  will  not  make  any  element 
nonzero  which  was  intentionally  set  to  zero  by  a  previously  ap- 
plied elementary  transformation.  Denoting  this  transformation  by 
Q  and  letting  TJ  be  the  completely  decomposed  matrix,  the  system 
is  solved  as  follows. 

Bx  =  y  (1.2.6) 

PBPPx  =  BPx  ■  y  (1.2.7) 


QcQpBPx  =  QcUPy  =  UPx  =  OcQpPx  (1.2.8) 


Px  =  «-1QcQpPy  (1.2.9) 
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x  «  PU"1QcOpPy  (1.2.10) 


Here  we  have  used  the  fact  that  PP  =  I  and  that  Pfc  *  P. 
The  transformation  U"1  is  performed  by  backsubstitution  on  U. 
Note  that  if  row  and  column  modifications  as  above  have  been  per- 
formed, the  algorithm  to  solve  the  system  does  not  change.  The 
relevant  multiplier  elements  in  the  elementary  transformations  of 
Qp  are  changed  only  if  row  modifications  were  performed,  while 
for  column  modifications  only  the  relevant  columns  of  0  have  dif- 
ferent  elements.   If  the  decomposition  is  terminated  so  that  the 

lower  right  hand  kxk  submatrix  remains  unfactorized,  the  modifi- 

1      2 

cation   of   i   rows   and  j  columns  requires  no  more  than  ^(i+j)n 

multiplications,  while  the  solution  of  the   new   system   requires 

2    11 
n   +  4k  .    If   the  factorization  of  the  kxk  submatrix  is  stored, 

the  subsequent  solutions  require  n  multiplications. 

To  give  the  reader  a  sense  on  how  large  k  can  be  such  that 

the  completion  of  the  factorization  for  each  new  solution  becomes 

comparatively  expensive  let  us  set  the  number  of  multiplications 

required  to  solve  the  system  equal  to  the  number  required  to  com- 

3 
plete  the  factorization,  that  is  n   =  ik   or  k=(3n)    For   n=100, 

k=45,  and  for  n=400,  k=113,  and  for  n=1000,  k=208. 
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CHAPTER  2 
A  SENSITIVITY  TRANSFORMATION  FOR  STUDYING  THE  PERTURBATIONS  OF 
SOLUTIONS  DUE  TO  THE  PERTURBATION  OF  THE  SYSTEM 
SECTION  2.1 

PERTURBATION  OF  A  SINGLE  ELEMENT 


The  first  example  to  be  investigated  is  the  change  of  the 
solution  due  to  the  change  of  one  element  in  the  system  of  equa- 
tions. This  result  was  first  deduced  by  Shermam  and  Morrison[8]. 
Though  the  result  here  is  the  same,  the  procedure  and  representa- 
tion are  different.  At  the  end  of  this  section  we  will  point  out 
how  this  is  so. 

Let  B  be  a  nxn  matrix  that  is  identical  to  the  matrix  B 
except  for  the  element  (i,j)  which  is  represented  as: 


Bij  ■  bij  +  6bij  •  (2.1.1) 


Let  B  be  decomposed  by  a  decomposition  technique.  Then  the  n 
columns  of  the  inverse  of  B  can  be  found  by  solving  the  decom- 
posed system  for  the  appropriate  unit  vector  as  a  right  hand 
side.   Suppose  that  x  is  the  solution  of 

Bx  =  y  (2.1.2) 

and  x  is  the  solution  of 

Bx  =  (B+6B)x  =  y  .  (2.1.3) 


22 


Premultiplying  the  above  equation  by  B   we  have 


(I-fB"16B)x  «  B-1y  *  x 


(2.1.4) 


where  6b  is  a  square  matrix  of  order  n  in  which  only  one  element 
is   nonzero,   that  is,  6bi-i»   Written  out  explicitly  (I  +  B-16B)  , 


is  of   the   form: 


Bii*bij 


Bj-i.i6bij 


l4»ji*ij 


bj*i.i6bij 


Bni6bij 


(2.1.5) 
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where  (B.  ,)  «  B  .   It  follows  directly  that 


kl 


x.  *  13 ,  (2.1.6) 

1  1+  5. .6b. 

ji  13 


and  for  Mj 


x.  »  xk  -  bki6b,  -x,  »  xk  -   — *i — iiJ-  .        (2.1.7) 
K    K    K1  1D  3    K     1  ♦  Bj^bij 


This  result  can  also  be  found   by   representing   the  per- 
turbed system  by 

(B  ♦  6b)x  *  (I  ♦  6BB-1)Bx  =  y  .  (2.1.8) 

If  we  denote  Bx*y  as  the  perturbed  system,  then  solving 

(I  +  6BB~1)z=y  (2.1.9) 

and  then 

Bx  =  z  (2.1.10) 

yields  the  solution  to  the  perturbed  system.   The  transforming 
matrix  (I  +  6BB~  '  is  of  the  form: 
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Csl 


to 


JO 


\o 


w 


la 


a 


in 
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-1 


Solving  (I  ♦  6BB~x)z  «  y,  for  Mi 


zk   yk  ,  (2.1.12) 


while   for    the   ifc     element, 


+    (1   ♦  *>ljBji)i1   +  6bi;jBji  +  1yi+1   +  ... 

...    +  <4bi;(Bjnyn    ,  (2.1.13) 


or 


yi  "  6bijBjiyi  +  6bijBj2y2  +  •••  +  6bijBji-iyi-i  + 

+    (1  +  ft^BJDii   +  »bi;jBji+1yi+1  +  ... 

...   +  «bi3Bjnyn    .  (2.1.14) 


Since 


n 


*1*,ijbjlt»*  =  6bijxj    '  (2'la5) 


(2.1.14)    reduces   to 
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yi   "  6bij\i   *    (1   *  6bijbji)zi   "  6bij5jiyi    •  (2.1.16) 


or 


(1  +  6biiBii)yi  -  ^iixi 
z. U_2i — i iJLJ.  (2.1.17) 

1  ♦  6bijBji 

6biixi 

■  y,  -    iJ-2 .  (2.1.18) 

1  +  6b.  .5. 


Therefore  z  can  be  written  as 


6biixi 

z  »  y  -    ( i-i-J )e.  (2.1.19) 

1  ♦  6bij5ji 

where  e.   is   the   unit  vector   of   the   i    component.    Since 
x  *  B~1z, 


*  -  x  - U-l B  1ei  .  (2.1.20) 


1  +  6b. .5 


Explicitly   for    Mj 


th.46b4.1X 


xk    *  x.     -       — !Si — LL_1_  (2.1.21) 


1   ♦  6bij5D] 


and   for    the   jtn   element, 
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Xj  * 


bii6bJiX3 
1  +  6bij5j. 


1  +  6bi-Bji 


(2.1.22) 


which  corresponds  to  (2.1.6)  and  (2.1.7). 

The  next  two  sections  will  consider  multiple  row  and 
column  perturbations.  We  will  see  that  the  analysis  of  such  per- 
turbations will  be  quite  straight  forward.  For  column  perturba- 
tions the  perturbed  system  of  equations  will  be  represented  as  in 
(2.1.4),  while  for  row  peturbations  the  representation  (2.1.8)  is 
appropriate.  In  Sherman  and  Morrison [8]  it  is  only  shown  that 
(2.1.6),  and  (2.1.7)  are  correct,  though  they  do  not  show  how 
they  arrive  at  the  results.  The  purpose  of  this  presentation  is 
to  demonstrate  a  procedure  which  can  be  used  to  derive  the  solu- 
tion for  an  arbitrary  perturbation,  rather  than  propose  a  solu- 
tion and  demonstrate  its  validity. 
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SECTION  2.2 

COLUMN  PERTURBATIONS 

In  this  section  the  method  just  described  is  extended  to 
columnwise  perturbations.  The  first  case  is  for  perturbations  in 
a  single  column  with  the  scaling  of  the  column,  while  the  second 
is  the  extension  to  two  columns. 

Without  loss  of  generality  let  us  suppose  that  the  first  k 
elements  of  the  j  column  are  to  be  perturbed  and  that  a  con- 
stant multiple  of  the  column  is  to  be  added  to  the  entire  column, 
that  is,  6b  is  an  nxn  matrix  which  has  all  zeroes  for  its  ele- 
ments  except  for  the  j   column,  where  the  j   column  is: 


J6blj+cbij  6b2j+cb2j   ...   *bkj+cbkj  cbk+lj   ...   cb^l*. 


(2.2.1) 


Again  the  matrix  (I  +  B~16B)  has  the  form  of  (2.1.5)   except  the 


j   column  is  replaced  by: 


where 


J**!      *2      '••      *j-l      <1+*j+c>      *j+l##'      ^nT  (2.2.2) 


'■  =  piiV^pj  ■  <2-2-3> 


and  (B..)*  B~  .   The  solution   is   found   oy   inspection  as   in 
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Chapter  2.1: 

and  for  i;*j# 

*i  "  xi  "  *ikj  (2.2.5) 

where  x  is  the  solution  of  (B  +  6B)x=y,  and  x  is  the  solution  of 
the  unperturbed  system  Bx*y. 

Now  suppose  that  two  columns  are  to  be  perturbed ,  say 
columns  j.  and  j2.  For  the  simplicity  of  demonstration  let  us 
suppose  that  only  the  first  k1  and  k2  elements  of  each  respective 
columns  are  perturbed.  So  6b  has  only  2  nonzero  columns,  and  the 
j1tn  is  of  the  form 

|6b1H   6b9_i    ...   6bk  •    0   ...    0  |fc         (2.2.6) 
I   1J1    ^1         KlJl  I 

and  the  j2  column  is: 

K  !t 

|6b,.    6b9.    ...   6bk  .    0   ...    0  r  .       (2.2.7) 


Now  two  columns  of  (I  +  B~16B)   are  not  unit  vectors,   namely 
columns  j.  and  j2.   The  jitn  column  will  be  denoted  by 


!*1      *2 


J1-1 


Dl 


j,-n#" 


*n| 


30 
(2.2.8) 


where 


*i   ■  pl^96bPh    ' 


(2.2.9) 


th 


and   the   j2       column  as, 


\#1     02      ••*      ^ 


jrl      (1+V      'j^l  —      0nj 


(2.2.10) 


where 


0.    *     2  5.    6b 
1        p-1   1P     P32 


(2.2.11) 


Again   by   inspection   the   j*        and   J2th   elements     of     x     are 


found  by  solving 


*|  1+0-i 


(2.2.12) 


Then    for    Mj1#    j2: 
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*k  *  xk  "  V^  "  *k8j2  •  (2.2.13) 


If  all  the  elements  of  a  column  are  to  be  perturbed  then 
it  is  necessary  to  have  the  entire  inverse  computed.  In  such  a 
case  it  would  be  wiser  to  use  the  factorization  modification  al- 
gorithm. In  many  cases  thouqh,  especially  in  Input  Output 
analysis,  only  a  few  elements  are  modified,  and  the  rest  are  just 
scaled.  In  such  cases  the  above  method  becomes  very  appealing 
since  the  simpler  closed  form  equations  provide  greater  insight 
into  the  system. 
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SECTION  2.3 

JKOW  PERTURBATIONS 

This  section  nearly  duplicates  the  presentation  of  the 
previous  section  except  that  the  perturbation  vectors  are  rows 
instead  of  columns. 

Without  loss  of  qenerality  let  us  suppose  that  the  first  k 
elements  of  the  j  row  are  to  be  perturbed  and  that  a  constant 
multiple  of  the  row  is  to  be  added  to  the  entire  row.  Then  6b  is 
an  nxn  matrix  which  has  all  zeroes  for  its  elements  except  for 
the  j    row.   The  j    row  is: 


6bjl*cbji   6bj2+cbj2   ••'   6bjk+cbjk   cbjk+l   •••   cbjn 


(2.3.1) 


i-l.  , 


Again  the  matrix  (I  ♦  6BB   )  is  of  the  same  form  as   (2.1.11)   of 
Chapter  2.1 ,   but  the  j    row  is  replaced  by: 


where 


I*!  *2      •"•   *j-l   (1+*j+c)   *j+l"-   *ni  '       (2.3.2) 


'"    =  pli6b3PV  '  (2-3'3) 


and  (15^)*  B"1.   By  inspection  the  solution   of   (  I  +6BB~1)z=  y 
is,  for  ij*j: 
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zi    ■   Yi  (2.3.4) 

while 


^V.  +    <*   +  *)    +  c»zj   ■  Vj  <2-3-5) 


^.-•A   +  ^Vj   -  ry,   +    (1   +  ^   +  c)zj  (2.3.6) 


pi^jp"?  "  Vj  +  (1  +  »*j  +  c,2j  ,2-3-7) 


and    therefore 


(J    +  >VyJ   •   D?/b3PXP 

zj  "    (1  * TC  H) ,2'3-8) 


where  x  is  the  solution  of  the  unperturbed  system   Bx=y.    If   we 
denote  z .  -  y^  as  t  then 


z  =  y  +  te.  (2.3.9) 


where  e.  is  the  unit  vector  of  the  j    component.    Denoting   the 
solution  of  the  perturbed  system  as  x ,    as  in  Chapter  2.1, 


34 
x  «  B"1z  *  B-1(y  ♦  te.)  *  x  +  tB-1e.  .  (2.3.10) 

Now  suppose  that  two  rows  are  to  be  perturbed,  say  rows  j, 
and  j2.  Again  purely  for  the  simplicity  of  demonstration  suppose 
that  only  the  first  k.  and  k2  elements  ot  each  respective  rows 
are  perturbed.  So  6b  has  only  2  nonzero  rows,  the  ji  of  which 
is  of  the  form: 


|6b.  ,   6b.  0      ...   6b.  .    0   ...    0  (2.3.11) 

I  *ll  31Z  D1K1  I 


and  the  j2  row  is: 


|6b.  ,   6b,  0      ...   6b,  „   0   ...    0  (2.3.12) 

32l  J2Z  32K2  I 


Now  two  rows  of  (  I  +  6BB   '  are  not  unit  vectors,  namely  rows  j, 
and  j2#   The  ji    row  will  be  denoted  by 

I  I 


1     ^2 


rl      (1+>V      ^j1+l-'-      *n 


(2.3.13) 


where 


ki 


*•»  =  pii'bi ,pBpi  '  ,2-3-14) 


th 


and  the  j2   row  as» 


where 


Sol' 


,-1 


ving  (  I  +6bb   ) z  ■  y,  for  i^j ,  j 


while  for  i=  Ji»Jo: 


where 


and 
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*1   *2   •'•   ^jrl   Cl^j  )   0jl+1 


•  •  •   0 


n 


(2.3.15) 


0,  -   2  6b.   T5  .  . 
1    p=l   D2p  pi 


(2.3.16) 


z  .  ■  y  . 
3        Y3 


(2.3.17) 


I  14*,     >*,    I  |z. 


I 


2   II  Jl| 

0^        1+0^      !z 


D2  M  3? 
I  I   2 


(2.3.18) 


y+     -  (l  ♦  A*  )y,  +  r*,  y-i  -  2  6b.  x 


3i  '3i     3o'3 


1   Jl     J2  J2    p=l    Jl 


JiP  P 


(2.3.19) 
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32 


Ya      -  (1  +  *a    )yi   +  0t  Vt   -   2   <*b.   xD  .        (2.3.20) 
J2         J2   D2     Dl  31    p=l    32p  p 

Redefining  t  to  be  z.   -  y.   and  defininq  u  to  be  z^   -  y+    ,  then 

31  31  D2    D2 

the  perturbed  solution  is: 


x  =  B"1z  =  B"1y  +  tB"1e.   +  uB*1e.   .  (2.3.21) 

Dl        J2 


If  an  entire  row  is  perturbed  then  only  its  associated 
column  of  the  inverse  is  needed  in  the  computation  of  the  per- 
turbed solution  vector.  The  computation  requires  2n+l  multipli- 
cations,  (column  perturbation  would  require  about  n  .) 

The  above  presentation  lends  particular  insight  into  the 
perturbation  of  row  elements  in  linear  equations.  When  several 
rows  are  perturbed  then  the  perturbation  of  the  solution  is  a 
linear  combination  of  the  associated  columns  of  the  nominal  in- 
verse. Therefore  we  can  make  a  rough  quess  on  whether  solutions 
are  sensitive  to  perturbations  of  several  rows  of  the  10  matrix 
by  inspecting  the  associated  columns  of  the  n^ninal  inverse. 
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CHAPTER  3 
A  NONLINEAR  INPUT  OUTPUT  MODEL 


As  we  have  seen  both  factorization  modification  and  solu- 
tion perturbation  have  comparative  advantages  in  the  solution  of 
modified  linear  equations.  Factorization  modification  has  the  ad- 
vantage of  simplicity  of  data  storage,  that  is,  when  a  row  or 
column  is  changed,  then  only  that  row  or  column  of  the  modified 
factorization   is   changed.    Also  the  repetitive  modification  of 

the  linear  equations  is  very  stable  numerically  since  the   refac- 

5 
torized  matrix   is  machine   identical   to  the  factorized  matrix 

when  the  modified  system  of  linear  equations  is  completely  decom- 
posed  from  scratch.   The  disadvantage  of  factorization  modifica- 

2  2 

tion  is  that  n  multiplications  must  be  performed,  and  n   matrix 

elements  must  be  read  in  from  secondary  storage.  In  some  cases 
solution  perturbation  requires  much  less  computation  than  factor- 
ization modification  to  achieve  the  same  result.  The  disadvan- 
tage of  solution  perturbation  is  that  if  many  elements  in  dif- 
ferent rows  are  to  be  changed,  as  in  the  case  of  modifying  all 
the  elements  in  a  single  column,  then  much  or  all  of  the  nominal 
inverse  must  be  computed.  Also  repetitively  computing  new  solu- 
tions due  to  changes  in  the  elements  of  the  matrix  is  difficult 
by  this  method.  In  this  section  we  will  utilize  the  comparative 
advantages  of  both  methods  in  solving  a  nonlinear  Input  Output 
model . 


5 

"Machine  identical"  means  that  the   bit   patterns   of 

the  two  matrices  are  identical. 
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There  has  been  very  little  work  presented  in  the  litera- 
ture about  nonlinear  10  models,  most  of  it  only  in  the  way  of  ex- 
istence proofs,  Sandberq[6)  and  Lahiri[?J,  and  V3ry  little  empir- 
ical application.  This  may  be  due  to  the  difficulty  of  obtaining 
data  on  the  nonlinear ity  of  the  transactions,  for  just  obtaining 
the  transactions  matrix  data  is  a  difficult  task  in  itself.  But 
the  nonlinear ity  of  the  transaction  elements  may  not  be  due  so 
much  to  the  nonlinearity  of  the  input  requirements  per  unit  out- 
put for  each  firm.  Rather  the  nonlinearity  may  be  due  to  the 
average  input  requirements  per  unit  output,  that  is  the  a^'s, 
shift  from  the  production  techniques  of  firms  that  cannot  in- 
crease output  towards  the  production  techniques  of  those  firms 
that  can.  A  good  example  of  this  type  of  nonlinearity  is  the 
transactions  to  the  oil  industry  if  oil  demand  would  change.  If 
the  demand  for  oil  by  consumers  and  industry  would  drop  drasti- 
cally then  the  transactions  for  drilling  equipment  and  real 
estate  would  drop  disporportionally  to  the  drop  in  total  oil  out- 
put since  oil  wells  would  last  for  a  longer  period  of  time. 
Given  more  time  there  would  be  fewer  chances  that  wells  that  are 
drilled  end  up  to  be  dry,  etc.  It  is  likely  that  at  lower  oil 
demand  more  oil  drilling  would  be  internally  financed,  and  thus 
the  amount  of   interest   charges  would  be  proportionally  less. 


The  transactions  matrix  is  AX  where  X  is  a 
diagonalization  of  the  total  demand  vector.  The 
transactions  matrix  is  the  interindustry  data  that 
governments  actually  collect,  the  A  matrix  is  then  derived 
from  the  transactions  matrix  after  the  total  demand  vector 
is  computed  by  summing  the  rows  of  the  transaction  matrix 
and  adding  the  final  demand  vector  to  the  resultant  vector. 
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Without  scarcity  driving  up  prices  firms  are  less  willing  to  take 
chances,  they  may  curtail  research  and  development  activity. 
Thus  the  composition  of  value  added  changes.  As  we  shall  see  in 
the  example  below,  if  the  demand  for  oil  would  rise  the  a. .  #s  for 
the  oil  industry  would  shift  toward  the  production  techniques  of 
the  exploratory  firms.  New  production  techniques  such  as  extrac- 
tion of  oil  from  shale  would  become  more  profitable,  and  average 
production  shifts  towards  this  technique.  By  classification  of 
firms  into  two  groups,  those  that  can  and  cannot  increase  output, 
we  have  an  elementary  procedure  by  which  we  can  represent  the 
nonlinearity  of  transactions. 

Let  j  be  the  index  corresponding   to  the  oil   industry. 
Suppose  that  oil  total  output  is  represented  by 

Xj  *  Xnj  +  X°j  '  (4.1.1) 


where  x . ,  xn • ,  x°^  are  total  demand  for  oil,  total  demand  for  oil 
supplied  by  new  oil,  and  total  demand  for  oil  supplied  by  old  oil 
respectively.  We  assume  that  only  the  producers  of  new  oil  can 
expand  the  output  of  oil  to  a  level  greater  than  x  ..  The  tran- 
saction from  the  i    industry  to  the  oil  industry  is 


ttj  =  aijXj  =  a°ijx°j  +  anij(Xj  -  x^)  .  (4.1.2) 


Therefore 
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/  o     .n   ,vo 
tii    <a  ii     ii,x  ii     n 
13*  r    x^  xj  r        ij         H.i.JJ 


for  x.^x  ^ .  (4.1.3)  exemplifies  the  shift  in  the  technical  coef- 
ficients mentioned  in  the  preceding  discussion. 

Though  oversimplified,  (but  the  linear  10  model  is  even 
more  simplified,)  this  example  does  illustrate  that  empirical 
work  in  this  area  may  be  more  practical  than  was  previously 
thought.  It  can  also  open  up  the  possibility  of  treating  invest- 
ment as  a  direct  input  to  production  instead  of  treating  it  as 
final  demand.  Investment  cannot  be  realistically  treated  as  a 
linear  function  of  total  output  since  if  existing  capital  is  used 
at  full  capacity  the  only  way  output  can  increase  is  by  invest- 
ment. Thus  investment  would  make  up  an  considerable  portion  of 
costs.  Alternatively  if  capital  is  not  being  used  at  full  capa- 
city then  there  will  be  little  investment  till  excess  capacity 
has  depreciated  away.  Observation  of  plant  capacity  levels,  and 
the  catagor ization  of  firms  in  a  industry  may  give  us  the  essen- 
tial information  needed  to  produce  a  viable  nonlinear  10  model. 

The  nonlinear  10  model  in  the  following  exposition  has  an 
restriction  in  the  form  of  the  nonlinear ities.  It  is  not  the  most 
general  model  that  could  be  presented,  though  the  formulation 
does  encompass  the  situations  that  we  would  most  likely  encounter 
in  the  real  world.  We  assume  that  the  nonlinearities  of  the 
tranasctions  are  only  a  tunction  of  the  buying  industries,  that 
is 
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tij  *  tiJ(XJ)  "  (4.1.4) 

This  is  a  crucial  assumption  in  this  algorithm  for  it  tremendous- 
ly reduces  the  order  of  the  computation.  Though  this  is  a  res- 
triction to  the  analysis,  it  is  difficult  to  visualize  an  example 
where  a  transaction  element  is  a  direct  function  of  total  demand 
other  than  the  buying  industry.  Even  in  the  literature  though 
this  assumption  is  seldom  made  or  utilized,  the  examples  present- 
ed usually  are  of  this  form.  Also  we  have  assumed  a  functional 
form  for  the  nonlinear  it ies.  Though  this  assumption  is  not  cru- 
cial to  the  exposition,  and  it  can  easily  be  relaxed,  we  shall 
see  that  this  representation  is  computationally  advantageous. 

A  truncated  power  series  will  be  used  as  the  functional 
form  of  the  nonlinearities.  The  use  of  only  three  terms  is  only 
for  the  simplicity  of  demonstration. 

x.  x. 

(4.1.5) 


where 


i  ■  «ij  +  Pij  +  'ij 


(4.1.6) 


Note  that  3i-j=ai-i  (x^)  .  We  denote  the  matrix  A"  and  vectors  x  and 
7  as  the  nominal  technical  coefficients  matrix,  »nd  nominal  total 
and  final  demand  vectors  respectively.   If  c<.  .=1,   and  Pi-\~Yi-hs^ 
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for  some  aj-i(x.i)»  then  the  function  is  constant,  as  in  the  linear 
10  model.  Varying  the  values  of  pi-  and  Yi-  from  zero  introduces 
nonlinear ity  in  the  technical  coefficient.  But  as  long  as  con- 
straint (4.1.6)  holds  then  the  nominal  output  vector  x  is  still 
the  solution  of  the  nonlinear  10  model  when  the  final  demand  vec- 
tor is  y. 

To  introduce  the  concept  of  the  model,  let  us  assume  that 
only  the  j  industry  has  nonconstant  technical  coefficients. 
These  technical  coefficients  are  only  a  function  of  the  total 
output  of  the  j  industry,  that  is,  x..  Let  y  be  a  vector  dif- 
ferent from  y  and  let  x1  be  the  solution  of 


(I  -  T^x1  =  y  .  (4.1.7) 


Clearly  there  is  little  chance  that  x1  is  the  solution   of 
the  nonlinear  equation 


(I  -  A(x))x  =  y  .  (4.1.8) 


Let 


AMx)  =  A(x)  -  A  .  (4.1.9) 


Note  that  AA(x)  is  nonzero  only  for  the  j    column.   By   solution 
perturbation 


(I  -  A(x))x  =  (I  -  A)(I  -  (I  -  AJ^AAfx))*  =  y     (4.1.10) 
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or 


(I  -  (I  -  A)_1^(x))x  =  x1  (4.1.11) 


where  x   is  the  solution  of  (4.1.7) 


The  jth  row  of  4.1.11  is 


n 


[1  -  I   5,-6  a.-fx,)]*,  =  x*  (4.1.12) 

i  =  l  J x    a j   J    J     J 


where 


-1 


(5^)  =  (I  -  A)"1  ,  (4.1.13) 


and 


(^aij(xj))  =AA(x)  .  (4.1.14) 


Since  the  x.  is  the  only  unknown  in  (4.1.12)  we  need  only  solve 
this  scalar  nonlinear  equation  for  x..  If  the  nonlinearities  are 
represented  as  a  truncated  power  series  as  in  (4.1.5) ,  then 
(4.1.12)  reduces  to 


X  •  X  • 

[1  -  en-*)  -  1)  -  a((-J-)2  -  l)]x   =  x* 
x.  I 
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(4.1.15) 


where 


n 


*   s  2   ^iiaiiPii  r  and  (4.1.16) 


e  =  2  B.x. y..  .  (4.1.17) 

i=l  D1  13   J 


Then  (4.1.15)  can   be   solved   by   a   numerical   method   such   as 
Newton's  method. 

The  advantaqe  of  the  functional  representation  (4.1.5)  is 
that  the  effects  of  the  nonlinear ities  of  all  the  elements  of  the 
j  column  of  A(x)  are  expressed  in  cr  and  0.  If  more  terms  are 
added  to  (4.1.5)  then  only  similar  terms  are  then  added  to 
(4.1.15).  Note  that  only  the  j  row  of  the  inverse  needs  to  be 
computed  to  solve  for  x . ,  but  once  we  have  solved  for  x .  we  now 
know  all  the  values  of  the  technical  coefficients  in  the  j 
column.  Using  factorization  modification  we  can  compute  the  en- 
tire solution  without  the  computing  the  inverse,  as  ve  would  have 
had  to  do  if  we  use  solution  perturbation  only.  During  the  fac- 
torization modification  we  can  simultaneously  compute  the  solu- 
tion x,  thus  saving  an  extra  pass  through  the  matrix. ^ 


3 

This  is  only  important   of   course   if   the   computer 

installation   does   not   have  enough  main  memory  to  hold  the 

entire  matrix. 
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The  extension  to  case  of  nonlinearities  in  more  than  one 
column  of  A  is  quite  straightforward.  If  k  columns  of  A  have  non- 
linearities,  then  the  k  corresponding  rows  of  (I  -  "R)~  must  be 
computed  and  a  k  order  nonlinear  matrix  equation  corresponding 
to  (4.1.15)  is  solved. 

This  algorithm  is  a  tremendous  savings  over  the  usual 
Newton's  method  solution  of  this  problem,  which  would  be 


xk+1  =  xk  -  (I  -  ^_T)"1[xk  -  A(xk)xk  -  y]  (4.1.18) 


k  t"  h 

where  x   is  the  approximate  solution  at  the  k    iteration,  and   T 
is  a  vector  valued  function  whose  j    element  is 


ti«  2   a..(x.)x.  .  (4.1.19) 


For     each     iteration     (4.1.18)     requires    approximately 

3        2 

n  /3  +  2n   +  ikn  multiplications  where  n  is  the  order,  i+1  is  the 

number  of  terms  in  the  truncated  power  series,  and  k  is  the 
number  of  nonlinear  columns.  The  method  proposed  here  requires 
only  about  k  /3  +  (i+l)k  multiplications.  If  n=100,  k=5,  and 
i=3,  then  an  iteration  of  (4.1.18)  requires  about  353,000  multi- 
plications, while  the  method  presented  here  requires  about  162 
multiplications.  Thus  the  proposed  method  has  a  computational 
savings  by  a  factor  of  2,180! 

Most  importantly  this  method  gives  us  a  very  simple  means 
to   compute  the  "marginal  inverse",  that  is  compute  the  incremen- 
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tal  total  output  due  to  an  increment  in  final  demand.  The  margi- 
nal inverse  is  the  inverse  of  the  Jacobian  matrix  evaluated  at  a 
solution.  Sandberq [6, Theorem  1]  has  shown  that  the  inverse  of 
the  Jacobian  evaluated  at  a  solution  will  approximate  the  pertur- 
bation of  solutions  around  the  solution  due  to  perturbations  of 
final  demand.  If  the  vector  c  is  the  perturbation  in  final 
demand,  and  the  vector  x^  is  the  actual  perturbation  of  total 
demand  then 

xp  =  J"1c  +  A»(c)  (4.1.20) 


where 


1  |Aj£j  !  !  ->  0  as   ||c||  ->     0  .  (4.1.21) 

The  norm  operator  is  the  Euclidean.   In  our  example  in  which  only 
column  j  is  nonlinear,  if  (Jik)  =  J,  then: 

jik   =   -aik,    for    k*i,j  (4.1.22) 

jRk    =    1    -   aRk^    for    Mj  (4.1.23) 

and 


Jij    =    "   aij    "    (xj)^5T7a(xj)    for    i^'  (4.1.24) 
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jjj  =  *  ~  *jj  "  (Xj)^T7a(Xj)  '  (4.1.25) 


Since  the  Jacobian  matrix  is  identical  to  the  10  matrix 
except  tor  the  j  column  we  can  use  factorization  modification 
to  replace  the  j  column  of  the  nominal  10  matrix  by  the  the  as- 
sociated column  of  the  Jacobian  matrix.  In  doing  so  we  have  ef- 
fectively computed  the  factorization  of  the  Jacobian. 

At  this  point  it  is  worthwhile  to  digress  to  a  fundamental 
result  of  10  analysis.  In  the  linear  10  model  price  is  computed 
by  the  identity 


(I  -  Afc)p  =  v   or   p  =  Afcp  +  v  .  (4.1.26) 


The  j    row  of  (4.1.26)  reads:  the  price  of  good   j,   that 
is,   p.   equals  the  average  unit  costs  of  inputs,  that  is,  (A  p) . 
plus  the  cost  of  value  added,  that  is,  v..   In  this  equation  pro- 
fit  rates  must  be  assumed,  and  the  costs  are  average,  not  margi- 
ns 

nal.  Assuming  nonincreasing  returns  suppose  that  we  are  given  a 
final  demand  vector  y  and  the  corresponding  total  demand  vector  x 
which  is  the  solution  of  the  nonlinear  10  model.  Then  we  can 
find  a  price  vector  p  such  that  it  is  profit  maximizing  for  each 
industry  to  produce  the  total  demand   vector   x   when   the   final 


o 

Concavity  of  the  technical  coefficients  functions  and 
value  added  functions  would  imply  nonincreasing  returns. 


4t 


aemand  vector  is  y.   Tne  equation  wnich  computes  tnis  price  level 
is : 


J  P  -   v 


(4.1.27) 


wnere  cue  i    element  ot  v"  is 


v(x.)  4,xi)5H-v(x1)  . 

1 


(4.1.2b) 


we  have  allowed  tor  nonlinear ities  in  value  added,  and  as  tor  the 
transaction  elements,  the  nonlinear ities  are  only  a  runction  ot 
that  particular  industry.  ol  course  tnis  definition  ot  value 
added  excludes  profits.  we  can  easily  verity  that  the  price 
vector  which  is  the  solution  ot  (4.1.27)  is  the  price  vector  that 
makes  x  the  profit  maximizing  output  vector.  tor  the  j 
industry  tne  per  unit  profit  is 


n 


2   a.  (x  •  )  -  v  (x- 
i=l   J   J       3 


4.1. 2b) 


'iherelore    total    profits    are: 


n 


"    3    "    VPJ    "    i!1aiJ(X3>     -    V(XD)]     • 


(4.1.3k,) 


xaKinq    the    derivative   with    respect    to    x      we    have 


j 


n 


cJx-77   j    s    ^   Pj    *    .21Piaij(xj)    +    v(x 


J     l-l  3  1 


)  J 


(4.1.31) 


or 


'j   '   J1Pi[aij(xj)    +    <xj>c£r*ij<*j>] 


V(Xj)      4      (Xj)^-V(Xj)      . 


(4.1.32) 


Notice  that  (4.1.24)  and  (4.1.2b)  are  the  general  form  ot  the 
elements  ot  the  Jacobian  matrix  J.  Also  the  right  hand  side  of 
(4.1.32)  corresponds  to  (4.1.2b).  by  differentiating  each  of  the 
industry's  profits  function  with  respect  to  that  industry's  total 
output  we  see  that  (4.1.26)  is  tne  first  order  condition  tor 
profit  maximization. 
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DISCUSSION  AND  CONCLUSIONS 

A  summary  ot  the  relative  advantages  and  disadvantages  of 
solving  modified  systems  of  equations  by  factorization  modifca- 
tion  using  elementary  transformations  and  by  solution  perturba- 
tion was  presented  in  the  beginning  ot  Chapter  3.  The  conclusion 
of  the  discussion  was  that  neither  method  had  a  clearcut  advan- 
tage over  the  other.  When  the  two  methods  are  used  in  conjunc- 
tion with  each  other,  as  in  the  nonlinear  10  model,  the  resulting 
algorithm  can  be  very  efficient. 

The  efficient  solution  of  10  equations  has  been  largely 
ignored  by  economists.  This  is  evident  by  the  fact  that  up  to 
now  only  large  batch  computer  systems  were  used  to  solve  10  equa- 
tions.  Using  the  algorithms  discussed  in  this  paper  a  368    ord- 

Q 

er  10  matrix  was  factor ized  on  a  minicomputer  installation.    The 
solution   computation  of   the   solution   for   an  arbitrary  final 

demand  vector  required  8  seconds  user  time  and  12  seconds  system 

nma  10 
time . 
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The  computer  was  a   Digital   Electronics   Corporation 

PDP11/50.  The  computer  was  operating  un^er  the  UNIX 
operating  system  developed  by  Bell  Laboratories.  The 
algorithms  were  coded  in  C,  the  principle  language  in  which 
much  of  the  UNIX  system  was  coded. 

10 

The  Unix  operating  system   indicates   two   types   of 

program   timinqs,   one   for   user   time,   which   is  the  time 

required  to  perform  the  actual   computation   in   the   user's 

program   area.   The  other  is  denoted  as  system  time  which  is 

the  computation  required  by  the  operating  system  to  perform 

principly   the   input  output  functions.   Since  a  368tn  order 

matrix  requires  approximately  1   megabyte   of   storage,   the 

computation  of  physical  block  addresses  of  the  data  consumed 

the  largest  portion  of  the  total  computation.    System   time 

would   be  significantly  reduced  by  performing  raw,  ie.,  pure 

direct  memory  access  (DMA)  input  output. 
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Double  precision  arithmetic  was  used  by  both  the  decompo- 
sition and  solution  programs.  As  a  test  of  the  numerical  stabil- 
ity of  the  algorithms  on  actual  data,  the  Bureau  of  Economic 
Analysis   (BEA)   368th  order  10  matrix  was  decomposed,  and  then  a 

system  of  equations  utilizing  this  matrix  was  solved.   The   resi- 

-14  .  . 

duals  were  of  the  order  of  10    when  double  precision  arithmetic 

was  used.   Thus  the  use  of   elementary   transformations   has   not 
been  detrimental  to  numerical  stability. 

As  a  test  of  the  nonlinear  10  model,  nonconstant  technical 
coefficients  were  introduced  into  five  columns  of  the  1967  368 
order  BEA  10  matrix.  The  iterative  solution  perturbation  algo- 
rithm required  3.3  seconds  user  and  6.1  seconds  system  time, 
while  the  factorization  algorithm  which  computes  the  entire  solu- 
tion vector  required  31  seconds  user  and  27  seconds  system. 
Since  the  solution  perturbation  algorithm  computes  the  total  out- 
put for  the  industries  which  have  nonlinear  input  technical 
coef f f icients ,  we  can  check  the  computation  of  the  factorization 
algorithm   by   comDaring   the   elements  in  the  solution  vector  to 

those  computed   by   the   solution   perturbation   algorithm.    The 

-14 
residuals  were  of  the  order  ot  10 

For  a  small  set  of  nonlinear  industries,  the  computation 
involved  in  the  solution  of  a  nonlinear  10  model  is  about  two  or 
three  times  the  computation  involved  in  solving  a  factorized  sys- 
tem of  equations.  In  terms  of  computation,  there  is  little  that 
bars   the   incorporation   of   the   algorithms   in   empirical    10 
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research.  It  may  be  futile  to  consider  the  construction  of  a  em- 
pirical nonlinear  10  model  with  all  ot  the  columns  being  non- 
linear, at  least  this  is  so  tor  the  present.  Many  times 
researchers  who  utilize  10  models  tocus  most  of  their  attention 
upon  one  or  two,  or  a  small  qroup  of  industries.  The  effects  of 
alternative  technical  coefficients  ot  a  small  group  of  industries 
are  often  analyzed.  ihe  framework  presented  could  easily  be 
molded  to  such  applications.  Finally  it  is  suggested  that  inter- 
industry data  collection  should  gear  some  ot  its  efforts  toward 
the  determination  ot  capacity  levels  of  the  individual  firms. 
Doing  so  will  allow  the  construction  of  the  nonlinear  investment 
functions  needed  to  make  the  endogenous  determination  of  invest- 
ment levels  in  10  models  possible. 
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